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Quantum-computing  ideas  are  applied  to  the  practical  and  ubiquitous  problem  of  fluid  dynamics  simulation. 
Hence,  this  paper  addresses  two  separate  areas  of  physics:  quantum  mechanics  and  fluid  dynamics  (or  specifi¬ 
cally,  the  computational  simulation  of  fluid  dynamics).  The  quantum  algorithm  is  called  a  quantum  lattice  gas. 
An  analytical  treatment  of  the  microscopic  quantum  lattice-gas  system  is  carried  out  to  predict  its  behavior  at 
the  mesoscopic  scale.  At  the  mesoscopic  scale,  a  lattice  Boltzmann  equation  with  a  nonlocal  collision  term  that 
depends  on  the  entire  system  wave  function,  governs  the  dynamical  system.  Numerical  results  obtained  from 
an  exact  simulation  of  a  one-dimensional  quantum  lattice  model  are  included  to  illustrate  the  formalism.  A 
symbolic  mathematical  method  is  used  to  implement  the  quantum  mechanical  model  on  a  conventional  work¬ 
station  The  numerical  simulation  indicates  that  classical  viscous  damping  is  not  present  in  the  one-dimensional 
quantum  lattice-gas  system. 
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I.  INTRODUCTION 
A.  Overview 

The  purpose  of  this  paper  is  to  show  that  a  phase-coherent 
quantum  computer  can  be  used  to  simulate  the  behavior  of  a 
system  of  massive  quantum  particles,  propagating  and  col¬ 
liding  on  a  discrete  space-time  lattice.  This  discrete  quantum 
particle  system  is  called  a  quantum  lattice  gas.  I  have  used 
principles  and  concepts  from  quantum  mechanics  instead  of 
from  classical  mechanics  to  formulate  “local  rules”  for  an 
artificial  microscopic  particle  dynamics.  In  a  quantum  lattice 
gas,  this  is  possible  because  a  network  of  two-energy-level 
quantum  systems  is  used  to  encode  the  configuration  of  par¬ 
ticle  occupancies  throughout  the  lattice. 

There  are  two  parts  to  this  paper.  First,  I  analyze  a  glo¬ 
bally  phase-coherent  and  entangled  quantum  lattice-gas  sys¬ 
tem  governed  by  the  many-body  Schrodinger  equation  of 
quantum  mechanics.1  The  many-body  Schrodinger  equation 
is  reformulated  as  a  Boltzmann  equation  of  kinetic  transport. 
Assuming  the  quantum  computer’s  wave  function  does  not 
decohere  by  uncontrolled  entanglement  with  the  external 
world,  the  main  analytical  result  of  this  paper  is  the  deriva¬ 
tion  of  a  lattice  Boltzmann  equation  that  exactly  describes 
kinetic  transport  at  the  mesoscopic  scale  in  the  quantum  lat¬ 
tice  gas.  That  is,  the  lattice-Boltzmann  equation  is  an  exact 
representation  of  the  particle  dynamics,  including  all  effects 
due  to  quantum  superposition  and  entanglement.  This  refor¬ 
mulation  of  many-body  quantum  mechanics  represents  a 
quantum  computing  application  geared  towards  the  direct 
simulation  of  physical  dynamical  models.  A  hydrodynamic 
fluid  simulation  is  considered  here  as  a  test  case. 

Second,  numerical  data  taken  from  an  exact  simulation  of 
a  globally  phase-coherent  quantum  lattice-gas  system  is  pre- 
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‘The  quantum  state  of  the  quantum  lattice  gas  is  said  to  be  glo¬ 
bally  entangled  when  qubits  in  the  system  are  entangled  with  other 
qubits  in  the  system  positioned  arbitrarily  far  away  in  the  lattice. 


sented.  The  simulation  method  uses  symbolic  mathematics  to 
implement  a  quantum  mechanical  system  in  the  second 
quantized  representation.  A  globally  phase-coherent  wave 
function  is  simulated  on  a  classical  computer.  This  is  pos¬ 
sible  because  the  number  of  spatial  sites  of  the  lattice  is 
small  and  the  number  of  qubits  pjer  site  is  few.  The  main 
finding  from  the  simulation  is  that  it  is  possible  for  mass- 
density  waves  to  oscillate  indefinitely.  The  simulation  con¬ 
firms  that  there  is  no  viscous  damping  in  the  hydrodynamic 
sound  mode  of  the  artificial  fluid,  i 

B.  Background 

Other  types  of  quantum  lattice  gases  have  been  studied, 
beginning  in  the  mid  1990s,  by  Bialynicki-Birula  [1],  Succi 
[2,3],  Meyer  [4,5],  and  Boghosian  and  Taylor  [6]  to  model 
the  relativistic  Dirac  equation  and  the  nonrelativisitic  Schro¬ 
dinger  equation.  In  contrast,  the  macroscopic  scale  behavior 
of  the  quantum  lattice  gas  presented  here  is  classical,  even 
though  the  microscopic  scale  dynamics  is  quantum  mechani¬ 
cal  rather  than  classical  in  nature,  lihe  quantum  lattice  gas 
reduces  to  a  classical  lattice  gas  only)  if  the  collision  process 
causes  a  particular  incoming  configuration  of  particles  to 
scatter  into  only  one  single  “outgoing”  configuration.2 

In  two  previous  papers  on  quantum  lattice  gases  [7,8],  I 
considered  a  quantum  spin  system  where  the  system  wave 
function  was  collapsed  into  a  tensor  product  state  over  the 
spins  (or  qubits)  after  each  collision  step.  This  allows  for 
local  entanglement  to  occur  temporarily  and  avoids  global 
entanglement  altogether  when  the  |  particles  propagate 
through  the  lattice  [7].  Allowing  for  only  short-range  and 
short-time  entanglement  of  qubits,  the!  quantum  lattice-gas 
system  is  described  at  the  mesoscopic  scale  by  a  lattice 
Boltzmann  equation,  with  a  local  collision  operator  that 
obeys  the  principle  of  detailed  balance  [8]  (we  may  refer  to 
this  model  as  a  factorized  quantum  lattice  gas).  It  provides  a 


2This  follows  since  it  is  a  direct  generalization  of  a  classical  lattice 
gas  with  quantum  bits  replacing  classical  bits: 
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way  to  implement  the  lattice  Boltzmann  equation  in  an  un¬ 
conditionally  stable  manner  on  a  classical  computer.  Al¬ 
though  quantum  mechanical  ideas  inspired  the  formulation 
of  the  collision  process,  in  the  end,  the  factorized  quantum 
lattice  gas  is  a  probabilistic  classical  process.  The  salient 
feature  of  the  factorized  quantum  lattice-gas  formulation  is 
that  it  is  suited  for  implementation  on  an  array  of  small  quan¬ 
tum  computers,  interconnected  by  a  classical  communication 
network.  Therefore,  the  previous  papers  do  not  address  the 
situation  where  quantum  superposition  and  entanglement  can 
spread  throughout  the  entire  quantum  computer.  This  situa¬ 
tion  is  treated  here. 

C.  Organization 

In  Sec.  II,  I  introduce  the  quantum  lattice-gas  formulation 
from  an  analytical  perspective.  The  quantum  lattice  gas  is 
treated  at  the  microscopic  and  mesoscopic  scales  in  Secs. 
II A  and  II B,  respectively.  When  the  quantum  computer  is 
fully  coherent  throughout  the  entire  course  of  the  simulation, 
the  collision  operator  is  nonlocal.  Evaluating  it  requires 
knowledge  of  the  entire  system  wave  function  on  the  quan¬ 
tum  computer.  An  exact  representation  of  the  quantum  lattice 
gas’  mesoscopic  behavior  is  developed  in  Sec.  II B.  Its  me¬ 
soscopic  behavior  is  governed  by  a  latticc-Boltzmann  equa¬ 
tion. 

The  quantum  lattice-gas  formalism  is  presented  from  a 
numerical  perspective  in  Sec.  III.  The  numerical  methodol¬ 
ogy  used  in  the  simulation  of  the  quantum  system  is  pre¬ 
sented  in  Sec.  Ill  A.  The  numerical  method  discussed  in  Sec. 
Ill  A  1  is  based  on  a  representation  of  a  universal  quantum 
gate  expressed  in  terms  of  the  creation  and  annihilation  op¬ 
erators.  The  symbolic  rules  used  to  carry  out  the  exact  simu¬ 
lation  is  described  in  Sec.  Ill  A  2.  A  simple  one-dimensional 
lattice-gas  model,  used  in  this  paper  for  test  purposes,  is 
described  in  Sec.  Ill  B.  I  have  included  various  computer 
simulations  with  both  classical  and  quantum  mechanical  mi¬ 
croscopic  dynamics.  The  classical  and  quantum  mechanical 
versions  of  this  simple  one-dimensional  lattice-gas  model, 
called  the  lD3Px  model,  are  described  in  Secs.  HI  B  1  and 
III  B  2,  respectively.  Simulation  results  are  presented  in  Sec. 
Ill  C.  The  classical  and  quantum  mechanical  simulations  re¬ 
sults  are  presented  in  Secs.  Ill  C  1  and  III  C  2,  respectively. 
The  classical  simulations,  provided  for  comparison  purposes, 
are  done  at  the  microscopic  scale  and  also  in  a  classical 
mesoscopic  mean-field  approximation.  Then,  I  present  an  ex¬ 
act  simulation  of  the  quantum  lD3Px  model,  with  three  qu¬ 
bits  per  site  for  small  systems.  Approximation  schemes  are 
needed  to  compute  the  many-body  dynamics  on  a  classical 
computer,  except  in  the  case  of  very  small  system  size  or 
systems  with  very  few  particles.  An  exact  quantum  simula¬ 
tion  of  a  small  cluster,  comprising  21  qubits,  is  carried  out  on 
a  conventional  workstation  using  a  symbolic  mathematics 
technique  that  is  described  in  Sec.  Ill  A.  The  numerical 
simulation  gives  us  a  way  to  understand  the  quantum  lattice- 
gas  method  in  concrete  terms  and  is  a  necessary  step  toward 
achieving  numerical  simulations  on  quantum  computers. 

A  brief  summary  of  the  results  and  a  few  closing  remarks 
are  given  in  Sec.  IV. 


II.  ANALYTICAL  TREATMENT 
A.  Microscopic  scale 

In  quantum  computing  [9,10],  a  two-level  quantum  bit 
(called  a  qubit )  represents  the  smallest  unit  of  information 
that  may  be  in  a  superposition  of  the  discrete  states  |0)  and 
jl).  A  qubit  |^)  =  a|0)+/8|l)  has  an  amplitude  a  of  being 
in  the  zero  state,  |0),  and  another  amplitude  /3  ofbeinginthe 
one  state ,  |1).  The  complex  coefficients  are  constrained  by 
|ar|2  +  j/?|2  =  l  so  that  the  probability  of  the  qubit  being  in 
the  zero  state  plus  the  probability  of  it  being  in  the  one  state 
is  unity.  For  any  unitary  quantum  computation,  one  can  de¬ 
scribe  the  algorithm  by  specifying  a  unitary  evolution  opera¬ 
tion,  in  our  case  formally  written  as  eiHrlh,  acting  on  the 
system  wave  function,  |''P(f))>  which  constitutes  the  state  of 
the  quantum  computer’s  “memory.”  With  N  qubits,  the 
quantum  state  |T,(r))  resides  in  a  large  Hilbert  space  with  2N 
dimensions.  A  new  quantum  state  |^(/  +  r))  is  generated  by 
application  of  a  unitary  operator  (which  could  be  represented 
by  a  unitary  matrix  of  size  2NX  2N)  for  a  short  duration  r  as 

|^(/+r))=e'"T/ft|^(r)>.  (2.1) 

By  repeated  application  of  eiHrlh,  an  ordered  sequence  of 
states  is  generated  and  each  one  is  (jiven  a  unique  time  label. 
If  the  first  state  is  labeled  by  t  then  the  next  one  is  labeled  by 
t+r,  and  the  next  by  t+ 2  r,  and  so  forth.  In  this  way,  think 
of  the  computational  time  advancing  incrementally  in  unit 
steps  of  duration  t.  Of  course  the  state  of  the  quantum  com¬ 
puter  exists  at  all  intermediate  times,  say  at  t+  t/2,  but  for 
our  purposes  we  need  to  consider  only  the  state  at  intervals 
of  the  time  step  r.  Formally,  the  quantum  computer’s  evolu¬ 
tion  is  invertible  by  application  of  the  adjoint  of  the  evolu¬ 
tion  operator 

\V{t-T))  =  e-iflTlh\V(t)).  (2.2) 

This  computational  picture  is  consistent  with  the  Heisenberg 
picture  of  quantum  mechanics.  For  any  reversible  algorithm 
chosen,  the  task  is  to  map  the  algorithm  onto  the  dynamical 
evolution  of  interacting  qubits  within  the  physical  device, 
which  can  be  driven  by  external  control. 

1.  Preliminaries 

Consider  a  quantum  computer  with  qubits  arranged  in  a 
lattice-based  array  with  the  following  properties: 

(1)  V  is  the  number  of  lattice  sites. 

(2)  B  is  the  number  of  qubits  per  site  (and  the  number  of 
nearest  neighbors). 

(3)  N—  VB  is  the  total  number  of  qubits. 

(4)  2N  is  the  size  of  the  full  Hilbert  space. 

(5)  2B  is  the  size  of  the  on-site  submanifold,  denoted  B 
(and  the  number  of  on-site  configurations). 
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TABLE  I.  Ket  symbols. 


Symbol 

Size  of  manifold 

Description 

m 

2n 

Total  system  ket 

!<a> 

2b 

On-site  ket 

k> 

2 

Qubit,  local  state  ket 

At  each  site  of  the  lattice  resides  a  group  of  qubits  acted 
upon  by  a  sequence  of  quantum  gates  [10-13],  whose  action 
is  mediated  by  external  control.  The  quantum  lattice  gas’ 
evolution  can  be  formally  expressed  as  a  special  case  of  Eq. 
(2.1)  where  e‘HTlh^SC  as  follows: 

|¥(*i,  . . .  ,xv;t+T))=SC I'T'Oti, . . .  ,xv;t)).  (2.3) 

In  Eq.  (2.3),  S  is  the  streaming  operator,  which  in  matrix 
representation  is  an  orthogonal  permutation  matrix  with 
components  being  either  0  or  1.  S  is  the  “classical”  lattice- 
gas  streaming  operator.  However,  in  Eq.  (2.3),  C  is  not  a 
classical  operator.  It  is  a  unitary  collision  operator.  In  gen¬ 
eral,  when  expressed  in  matrix  form,  C  has  complex  compo¬ 
nents.  (The  quantum  lattice  gas  reduces  to  a  deterministic 
classical  lattice  gas  if  t  is  a  permutation  matrix  with  0  or  1 
components.  If  and  when  C  is  stochastically  switched  be¬ 
tween  different  permutation  matrices  during  the  dynamical 
evolution,  then  the  quantum  lattice  gas  reduces  to  a  probabi¬ 
listic  classical  lattice  gas.)  Finally,  in  Eq.  (2.3),  I  have  ex¬ 
plicitly  labeled  the  wave  function’s  dependence  on  all  the 
coordinates  of  the  lattice  to  emphasize  that  the  wave  function 
is  a  lattice-based  field  quantity. 

In  general,  the  operator  C  can  cause  mixing  of  outgoing 
collisional  configurations  at  each  site  of  the  lattice,  locally 
entangling  the  qubit  states  within  a  lattice  cell  of  size  /.  The 
operator  S  then  causes  particles  to  move  from  one  site  to  the 
next,  by  exchanging  qubit  states  between  nearest  neighbor¬ 
ing  sites.  Although  the  application  of  S  causes  the  particles 
to  move  just  as  they  would  in  the  streaming  phase  of  a  clas¬ 
sical  lattice  gas,  it  also  causes  global  superposition  and  en¬ 
tanglement  of  all  the  qubit  states,  if  local  entanglement  has 
already  been  caused  by  C.  In  this  way,  quantum  entangle¬ 
ments  are  spread  throughout  the  lattice  by  the  action  of  §. 

I  will  use  the  following  convention  for  indices. 

(1)  Small  roman  letters  ( a,b,c )  for  the  momentum  direc¬ 
tions  on  the  lattice,  a  e  {0, . . .  ,B  - 1 }. 

(2)  Greek  letters  (a,p,y)  for  specifying  qubits,  a 

e{0 . N- 1}. 

(3)  Middle  roman  letters  (ij,k)  for  the  spatial  dimen¬ 
sions,  *e{l, . . .  ,D}. 

2.  System  wave  function 

Let  |T'),  \>p),  and  \q)  denote  the  total  system  ket,  on-site 
ket,  and  qubit  ket,  repectively,  as  shown  in  Table  I.  The 
quantum  computer’s  total  wave  function  can  in  general  be 
expressed  as  a  linear  combination  of  tensor  product  states 
over  all  the  qubits 
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|^(x,,  ...,xv\t))=  2  A(<h . 4n) 

ill 

X|?i)®  •  •  -®kw)-  (2.4) 

Here  the  summation  indices  qa  are  either  zero  or  one,  for 
l^a=sA.  Each  tensor  product  \q\)®  ■  •  is  a  basis 

state  and  |1F)  is  a  pure  classical  state.  The  number  represen¬ 
tation  (2.4)  is  used  in  the  numerical  quantum  simulation  pre¬ 
sented  in  Sec.  Ill  C.  I  would  like  to  establish  a  convention 
for  representing  the  system  ket  as  a  linear  combination  of 
tensor  product  states  that  are  lattice-site  specific.  Let  \tp) 
denote  an  on-site  ket  formed  over  the  qubits  at  a  single  site 
of  the  lattice 

I  <Kx,t))=  X  a(qx,,..,qB) 

ki . 

Xki(x,/))®- •  -®\qB(x,t)).  (2.5) 

The  system  wave  function  (2.4)  can  in  general  also  be  ex¬ 
pressed  as  a  linear  combination  of  tensor  product  states  over 
all  the  on-site  kets 

|^(xi . xv;t))=  X  A(i/ri , . . .  ,ifry) 

\<h . (M 

X|^i)®- (2.6) 

where  the  shorthand  notation  \4>„)^\4'(xn  ,t))  is  used.  Here 
the  indices  ip„  (for  1  K)  in  the  sum  represent  the  num¬ 
bered  basis  states  in  the  on-site  manifold  B.  So  they  are  in 
the  range  0*s  ipn^2B~l .  The  coefficients  A  account  for  all 
the  global  superpositions  between  lattice  sites. 

3.  Unitary  collision  matrix 

Collisions  are  implemented  independently  at  each  site  of 
the  lattice.  Hence,  all  sites  can  be  collided  in  parallel,  homo¬ 
geneously  across  the  entire  system.  The  collision  operator  C 
is  therefore  expressible  in  tensor  product  form  since  local 
quantum  superposition  of  outgoing  on-site  configurations  oc¬ 
curs  only  within  each  25-dimensional  submanifold  B.  The 
2nX  2n  collision  matrix  t  can  be  written  as  the  following 
tensor  product: 

v 

£=  ®  £>,  (2.7) 

*=i 

where  the  on-site  collision  matrix  U  is  a  2BX2B  unitary 
matrix.  It  acts  on  the  on-site  ket 

\<fi'(x,t))  =  U\if/(x,t)).  (2.8) 

The  prime  on  the  left-hand  side  (LHS)  of  Eq.  (2.8)  indicates 
that  the  ket  is  an  outgoing  collisional  state.  Using  the  repre¬ 
sentation  (2.6)  of  the  system  ket,  the  postcollision  system  ket 
is 
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!¥'(*!, . . .  ,xv;t))  =  C .  ..,xv;t)) 

X|^j)0  •  -  ■®U\lpy) 

Wl . <p'y) 

X|^)®--.®|  tf>'y).  (2.9) 

An  equivalence  class  is  defined  as  a  set  of  basis  states  that 
correspond  to  particle  configurations  with  the  same  mass  and 
momentum  (and  energy  if  that  is  also  defined  in  the  lattice- 
gas  model).  The  on-site  unitary  collision  operator  0  acting 
on  the  5-submanifold  itself  is  block  diagonal  over  the 
equivalence  classes.  Consider,  for  example,  the  quantum 
lD3Px  lattice  gas  (see  Sec.  HI  B  1  for  a  detailed  description 
of  the  lD3Px  lattice-gas  model).  There  are  two  conserved 
quantities  for  this  one-dimensional  system:  the  mass  and  the 
momentum  along  the  x  axis.  Hence,  there  is  only  one  equiva¬ 
lence  class  and  it  has  two  members,  a  two-body  head-on 
configuration  and  a  configuration  with  a  single  rest  particle. 
Both  configurations  have  m= 2  and  p— 0.  The  equivalence 
class  is  comprised  of  the  following  on-site  kets: 

|3)  =  |011>, 

|4)  =  |100>. 

A  general  outgoing  ket  in  this  mass-momentum  sector  of  the 
on-site  submanifold  is  a  linear  combination  of  these  two, 

a|01 1)  +  /?|  100),  (2.10) 

where  a  and  5  are  complex  numbers.  So  the  collision  matrix 
0  for  this  one-dimensional  quantum  lattice  gas  has  one 
block.  It  has  a  2  X  2  block  for  mixing  the  head-on  and  rest 
particle  configurations.  In  general,  0  is  block  diagonal  over 
the  equivalence  classes  [7],  Each  block  of  0,  associated  with 
an  equivalence  class  of  size  n,  is  a  member  of  the  U(n) 
unitary  group. 

B.  Mesoscopic  scale 

1.  Occupancy  probability  and  the  mass  and  momentum  densities 

The  probability  of  occupancy  at  time  t  of  the  ath  local 
state  is  denoted  fjt).  Let  the  ath  local  state  be  associated 
with  a  displacement  vector  ea  at  position  x.  Also,  let  na 
denote  the  number  operator  for  the  ath  local  state.  That  is, 
«0|^(/))  has  eigenvalue  1  or  0  corresponding  to  the  ath 
local  state  being  occupied  or  empty  at  time  t.  A  fundamental 
construct  of  the  quantum  lattice-gas  formalism  is  that  the 
probability  of  occupancy  fa(t)  is  expressed  in  terms  of  the 
quantum  mechanical  density  matrix  6(/)=|>lr(f))0P(0l  as 
the  following  trace: 


fa(t)=fa(x,t)=Tr[e(t)nal  (2.11) 

In  the  literature  on  classical  lattice  gases  and  the  lattice- 
Boltzmann  equation,  fa(x,t)  is  referred  to  as  the  single¬ 
particle  distribution  Junction ,  and  it  is  defined  at  the  mesos¬ 
copic  scale.  For  classical  lattice  gases,  numerical  estimates  of 
fa(x,t )  are  obtained  either  by  ensemble  averaging  over  many 
independent  microscopic  systems  or  by  coarse-grain  averag¬ 
ing  over  space-time  blocks  with  a  single  microscopic  system. 
For  the  quantum  lattice  gas,  the  fa(x,t)  is  the  expectation 
value  of  the  operator  na  determined  by  repeated  measure¬ 
ment  of  single  microscopic  realizations  or  by  direct  measure¬ 
ment  of  an  ensemble,  as  occurs  in  nuclear  magnetic  reso¬ 
nance  quantum  computers  [14,15],  So  the  definition  (2.11) 
also  defines  fa(x,t )  at  the  mesoscopic  scale. 

Let  aQ  denote  the  first  local  state  within  the  group  of 
local  states  at  position*  of  the  Bravais  lattice.  In  addition,  let 
aQ  correspond  to  the  displacement  vector  e0.  Next,  suppose 
the  local  states  are  numbered  in  a  systematic  and  well- 
ordered  fashion  so  that  each  local  state  a=ao+a,  for  all 
a  e  {0,1, . . .  ,B—  1 },  resides  at  position  x.  Note  that  with  this 
numbering  scheme,  the  directional  index  a,  associated  with 
the  ath  local  state,  is  found  by  the  modulus  operation  a 
=  (a  modi?).  Then,  the  local  mass  density  and  the  momen¬ 
tum  density  at  *  and  t  can  be  expressed  in  terms  of  the 
occupancy  probability  fa(x,t)  following  the  convention  used 
to  define  the  mass  and  momentum  densities  in  a  classical 
lattice  gas 

B  aO+B 

p(x,t)=  lim  'Z  mfa(x,t)  =  lim  2  wTr[e(r)nJ, 


B 

p(x,t)Vi(x,t)=  lim  Z  mc2eafa(x,t) 

/,-o 

aQ+S 

=  lim  2  rnc2e{amodB)iTr[e(t)nal 
*o  “=“o 

(2.13) 

The  mass  and  momentum  densities  are  considered  “macro¬ 
scopic”  field  quantities.  They  are  only  well  defined  in  the 
continuum  limit,  where  the  primitive  cell  size  of  the  lattice 
approaches  zero.  However,  for  practical  considerations,  they 
are  approximated  by  high  resolution  grids  with  small  but 
finite  cell  size. 

To  experimently  determine  the  mass  density  or  momen¬ 
tum  density  at  a  site  *  at  time  t  in  an  actual  quantum  system, 
it  is  necessary  to  know  the  probability  of  occupancy  of  all 
the  local  states  at  that  site  fa(x,t)  for  a—  I, . . .  ,B,  according 
to  the  definitions  (2.12)  and  (2.13).  The  probability  of  occu¬ 
pancy  fa(x,t)  of  each  local  state  depends  on  the  polarization 
of  the  corresponding  qubit  |  qa(x>t))-  aa(x,r)|0) 
+5fl(*,0|l).  A  Von  Neuman  measurement  of  the  state  of 
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Qubits 

k> 

W) 

Local  state 

a 

af 

Position 

X 

x'  =x+/ea 

Momentum 

p  =  mcea 

P'-P 

this  qubit  will  yield  a  value  of  either  0  or  1,  with  probability 
\aa(x,t)\2  or \f3a(x,t)\1,  respectively,  since  the  measurement 
causes  a  collapse  of  the  quantum  wave  function.  A  single 
value  obtained  by  this  stochastic  measurement  process  is  not 
sufficient  to  determine  fa(x,t).  Therefore,  to  obtain  an  esti¬ 
mate  of  the  expected  equilibrium  values  of  the  mass  and 
momentum  densities,  it  is  necessary  to  either  ensemble  aver¬ 
age  over  many  realizations  of  the  microscopic  system  or 
coarse-grain  average  over  space-time  blocks  within  a  single 
microscopic  realization.  In  this  regard,  the  amount  of  effort 
needed  to  obtain  estimates  of  the  densities  is  identical  for  the 
quantum  system  and  classical  lattice-gas  sysytems.  A  quan¬ 
tum  computer  that  provides  a  direct  means  for  measuring  the 
expected  state  of  a  qubit  (such  as  is  possible  with  an  NMR 
quantum  computer)  would  be  a  more  natural  choice  for 
implementing  this  quantum  lattice-gas  algorithm. 

If  measurements  were  made  at  each  and  every  site,  and  at 
every  time  step  of  the  dynamics,  then  the  quantum  lattice-gas 
system  is  effectively  “factorized”  in  such  a  way  that  the 
quantum  computer’s  wave  function  is  always  collapsed  into 
a  tensor  product  state.  This  type  of  factorized  quantum 
lattice-gas  simulation,  with  continual  and  homogeneous  mea¬ 
surement  of  the  qubits,  results  in  a  probabilistic  classical 
lattice-gas  simulation  [8].  Yet,  even  in  this  case,  the  value  of 
the  transport  coefficients  can  differ  from  those  of  the  classi¬ 
cal  lattice  gas. 


2.  Mesoscopic  transport  equation 
Let  us  consider  two  qubits  | q)  and  | q'),  which  are  located 
at  neighboring  sites  x  and  x'  -x+/fea ,  respectively.  I  shall 
refer  to  the  local  states  encoded  by  these  two  neighboring 
qubits  by  their  numerical  labels  a  and  a',  respectively. 
Next,  suppose  these  local  states  may  be  occupied  by  particles 
with  momentum  mcea .  Following  this  construction,  the  ac¬ 
tion  of  the  streaming  operator  S  causes  a  particle  to  move 
from  site  x  to  the  neighboring  site  x' ,  hopping  from  local 
state  a  with  momentum  p=mcea  to  local  state  a'  with  the 
same  momentum  p '  -p.  This  labeling  convention  is  summa¬ 
rized  in  Table  II.  With  this  understanding,  we  can  write  the 
identity. 


(2.14) 

This  is  a  simple  mathematical  way  of  stating  the  following: 
If  you  make  a  measurement  of  the  occupancy  of  local  state  a 
before  streaming,  the  result  you  get  must  be  the  same  as 
when  you  make  a  measurement  of  a'  after  streaming. 

The  first  step  toward  deriving  a  microscopic  transport 
equation  for  the  quantum  lattice  gas  is  to  rewrite  Eq.  (2.3)  as 


(^(/)|^^t|^(f+T))  =  (^(,)|^ac|^(t)), 

(2.15) 

which  is  done  by  multiplying  through  from  the  left  by 
(^(/)|C^/!a51',  and  then  using  the  identity  5^5=1.  From  the 
identity  (2.14),  we  know  that  nJP—f$na< .  Using  this  fact 
in  the  above  equation  allows  us  to  write  it  as 

(^(t)|Ct5t«£t,|^(t+r))  =  (^(0|Ct«aC|^(/)). 

(2.16) 

The  “bra”  vector  on  the  LHS  of  this  equation  can  be  sim¬ 
plified  using  the  adjoint  of  Eq.  (2.3),  which  is  (^(/+t)| 
=  (^F(0|CtSt,  so  that  we  obtain  the  following  result: 

(^(f+T)|na)|^(/+r))  =  ('^(r)|C't«aC|^(0). 

(2.17) 

Using  Eq.  (2.11)  and  referring  to  Table  II,  Eq.  (2.17)  ex¬ 
presses  the  probability  of  occupancy  of  local  state  a’  at  site 
x+fsea  at  time  t+  r  in  terms  of  a  matrix  element  evaluated 
at  the  neighboring  site  x  and  at  the  earlier  time  t.  That  is, 

/a(;+/vea,t+r)  =  <^(t)|CtnoC|^(0).  (2.18) 

We  may  add /fl(x,f)-('SF(0l««|^r(0)»  which  vanishes  by 
definition,  to  the  right-hand  side  (RHS)  of  Eq.  (2.18).  Then, 
we  recognize  Eq.  (2.18)  is  a  transport  equation  for  the  par¬ 
ticle  occupancies.  The  result  is  a  lattice-Boltzmann  equation 
for  the  quantum  lattice-gas  system 

fa(x+/Ja  ,t+  T)^fa(x,t)  +  nrw,  (2.19) 

where  the  collision  term  is  expressed  as  the  following  matrix 
element: 


n™so(V)  =  {V(t)\frnaC-na\y(t)).  (2.20) 

An  alternative  derivation  of  Eq.  (2.20),  carried  out  in  the 
continuum  limit,  is  given  in  Appendix  A.  In  practice,  we  will 
not  be  able  to  analytically  evaluate  Eq.  (2.20)  for  a  large 
quantum  lattice-gas  system  with  global  entanglement  be¬ 
cause  of  the  exponential  size  of  the  |^)  ket.  However,  it  is 
possible  to  formally  express  the  collision  term  0™eso  when 
1^)  is  represented  as  the  linear  combination  (2.6).  This  is 
done  as  follows: 


a™es0= 


{^i 


2  2  , 

. M  {^1 . 


.  .  ,ipy) 


XA(ipiy  ...,t/rv)(t/t i-|®- 


®(<M<?tna6-ndv!'i)®  •  •  (2.21) 

Moreover,  it  is  possible  to  express  Q™eso  in  terms  of  the 
on-site  number  operator  na,  which  is  represented  by  a 
2bX2b  matrix.  That  is,  na  acts  only  in  the  submanifold  B  on 
the  qubits  at  a  single  site.  We  write  the  A-qubit  number 
operator  na  as  a  K-fold  tensor  product  that  has  a  single 
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5-qubit  number  operator  na  located  at  the  nth  site  index 
corresponding  to  the  position  vector  x„  as 

•  •  -®1,  (2.22) 

where  1  denotes  the  2BX2B  identity  matrix.  The  collision 
operator  C^naC—na ,  can  then  be  written  as 

\®\®---®{U^naU-na)®---®l=®vx=lCla, 

(2.23) 

where 


U'naU-na, 

1, 


otherwise. 


(2.24) 


Using  Eqs.  (2.7),  (2.22),  and  the  orthonormality  of  the 
on-site  kets  (if/„'\if/n)  =  8„,„,  Eq.  (2.21)  reduces  to  a  local 
matrix  element  evaluated  at  single  site  xn<=xn—x  as 


Then  using  Eq.  (2.22),  the  lattice-BoItzmann  equation  for  the 
quantum  lattice-gas  system  becomes  a  local  equation  that 
can  be  easily  simulated  on  a  classical  computer  [7,8]. 

3.  The  approach  to  steady-state  equilibrium 

The  system  is  said  to  be  in  steady-state  equilibrium 
(which  may  also  be  called  thermodynamic  equilibrium )  when 
the  system  ket  |'P'eq(r))  is  an  eigenvector,  with  unity  eigen¬ 
value,  of  the  collision  operator  C, 

Cl^H^eq).  (2.30) 

The  value  of  the  probability  of  occupancy  (2.1 1)  is  then  de¬ 
termined  from  |’Peq)  as 

/^)  =  (^(r)|«J^(r)>.  (2.31) 

Notice  by  the  definition  (2.30)  for  steady-state  equilibrium, 
the  collision  term  (2.20)  in  the  lattice-Boltzmann  equation 
vanishes, 


ar°=2  2 

K'  {0i.---.0r} 

XA*(lpl . +  U  ...,lpy) 

. */'v)(fa'\OinaO-na\iftH). 

(2.25) 


Let  us  make  the  following  definition: 


K(>Pn' M- 

{01 


. . ,<Pn  + 1  >  •  •  •  ,<Pv) 

XAifa . >//„ . tfty).  (2.26) 


The  quantity  ,<p„)  represents  the  superposition  of  the 
on-site  basis  states  at  site  x  with  all  the  other  on-site  basis 
states  in  the  system  at  the  other  sites.  With  this  definition, 
Eq.  (2.25)  can  be  written  in  a  simpler  way. 


ftr°=2  2  Wn'M(<Pn'\0^aU-nM. 

0»-  0» 

(2.27) 

If  each  on-site  state  is  not  entangled  or  superposed  with  any 
other  on-site  state,  then  TZ  can  be  written  in  factorized  form, 
TZ(  tp„ /  ,ip„)  =  C(4/„')C(<J/„)-  In  this  case,  Eq.  (2.27)  is  simpli¬ 
fied, 


nr°=(>P\0^aU-na\>p),  (2.28) 

where  the  coefficients  C(ipn)  specify  any  local  superposition 
and  entanglement 


C(<pn)\ >!>„)■  (2.29) 

0/1 


fln,es°(|^eq))  =  0  (2.32) 

Therefore,  at  steady-state  equilibrium,  the  occupancy  prob¬ 
abilities  are  unchanging  over  time.  That  is,  |1Feq)  is  the 
ground  state  of  the  system.  The  distribution  along  the  mo¬ 
mentum  directions  of  the  particle  occupancies  are  uniform, 
so  the  local  configurations  are  perfectly  symmetric,  and 
0“eso  cannot  cause  any  further  changes. 

HI.  NUMERICAL  TREATMENT 
A.  Methodology 
1.  Universal  two-qubit  gate 

In  this  section,  I  write  a  two-qubit  universal  gate  in  terms 
of  the  creation  and  annihilation  operators  of  the  second  quan¬ 
tized  formulation  of  quantum  mechanics.  A  classic  text  on 
second  quantization  is  by  Fetter  and  Walecka  [16].  For  our 
purposes,  the  two-qubit  gate  is  a  member  of  the  special  uni¬ 
tary  group  SU(2);  I  neglect  the  overall  phase  factor  because 
this  does  not  affect  the  quantum  lattice-gas  dynamics.  If  0  is 
a  member  of  SU(2),  it  can  be  parametrized  using  three  real 
numbers  £,  f,  and  0  as  follows: 

„  /  e'^cosd  -c,fsin0\ 

^  \—  e~‘lsin0  ~e~'^cos  0/’ 

Let  a\  and  att  denote  the  creation  and  annihilation  operators 
for  ath  spin  of  a  fermionic  quantum  spin  system.  Then  the 
spin- 1  creation  and  annihilation  operators  satisfy  the  anti¬ 
commutation  relations 

{aa,al}=8ap,  (3.2) 

{aa,ap}  =  0, 

{al,al}  =  0. 
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The  spin  number  operator  n^a\aa  has  eigenvalues  of  1  and 
0  in  the  number  representation  when  acting  on  a  pure  state, 
correspondirig  to  the  ith  spin  being  up  and  down 

sz=  -  j,  respectively. 

Consider  a  fermionic  spin  system  with  a  total  of  N  spins 
whose  system  ket  is  denoted  by  |M').  Acting  on  this  system 
ket  with  a  unitary  operator,  we  would  like  to  entangle  the 
two  spin  states,  the  states  of  the  a-th  and  /3th  spins,  according 
to  the  components  of  the  special  unitary  matrix  (3.1).  Let 
Yap  denote  a  square  2NX  2N  matrix  that  does  this.  I  define 
Yap  in  terms  of  the  multispin  creation  and  annihilation  op¬ 
erators  as  follows: 


where  X  is  a  symbol  used  to  denote  what  I  call  the  null  state 
that  accommodates  Pauli  exclusion  and  destruction  on  the 
vacuum.  That  is,  the  symbols  and  a  represent  the  single¬ 
spin  (or  single  qubit)  creation  and  annihilation  operators,  re¬ 
spectively. 

Next,  all  the  basis  states,  in  the  number  representation,  are 
encoded  by  the  symbol  ’Ff.s],  where  0^ssS2N—  1,  for  a 
system  with  N  spins.  That  is,  the  states  are  binary  encoded 
and  labeled  by  A-bit  integers.  The  state  1P[0]  is  called  the 
vacuum  state.  The  symbolic  rules  embodying  the  multiple- 
spin  creation  and  annihilation  operators  are  defined  in  terms 
of  the  single-spin  rules 


Y0j8=l+e-'^sin  6  a^a  a+ e‘hm  8a\a^-{  1  Telcos  8)na 
-(1  -e~'^cos  0)np—2i  sin£cos  8nanp  (3.3) 
for  a  A  j8.  Its  matrix  representation  for  a  two-qubit  system  is 


I1 

0 

0 

° 

■  0 

e'^cos  8 

— e'^sin  8 

0 

Y= 

0 

-e~lism  0 

~£-,fCOS  0 

0 

[ 0 

0 

0 

-1/ 

In  Appendix  B,  I  demonstrate  why  Yafl  is  manifestly  unitary 
and  an  appropriate  formulation  of  a  universal  quantum  gate. 

In  the  special  case  when  8-tt/2,  i=0,  and  £=0,  then 
Y ap  reduces  an  interchange  operator 


(-l)-s«at[(5A2“)=s-a] 

.0, 


if  j=X, 

(3.12) 


=X, 
(3.13) 

where  0*Saf=sA— 1  and  where  the  factor  (— l)s“  appearing 
in  Eqs.  (3.12)  and  (3.13)  accounts  for  a  phase  change  of  rr 
radians  induced  by  commuting  spins.  In  the  number 
representation  each  basis  state  is  denoted  by  a  ket 
\nln2-  ■  -na-  ■  -nN),  where  each  n  is  either  1  (particle 
present)  or  0  (no  particle  present).  The  phase  factor  Sa  is 
defined  by 


a[a-, ¥[$]]= 


(-1)^u[(sA2")=*«] 

0,  if 


Xap-l  +  al<ta+aiap-na-np,  (3.5) 

which  is  a  NOT  gate  (see  Appendix  B). 

2.  Symbolic  mathematics  method 

It  is  possible  to  simulate  the  exact  quantum  mechanical 
evolution  of  a  quantum  spin  system  using  computational 
symbolic  mathematics.3  To  test  the  quantum  lattice-gas 
method,  I  implemented  the  algorithm  using  version  4  of 
MATHEMATICA  [17].  Letting  1  and  0  represent  spin  up  and 
spin  down,  respectively,  the  first  step  is  to  define  a  set  of 
rules  that  encode  the  Fermionic  anticommutation  relations 
(3.2) 


at[0]=l, 

(3.6) 

at[l]=X, 

(3.7) 

at[X]=X, 

(3.8) 

a[0]=X, 

(3.9) 

a[  1]=0, 

(3.10) 

a[X]  =  X, 

(3.11) 

3I  developed  this  symbolic  method  in  1991  at  Brandeis  Univer¬ 
sity,  see  http://xyz.plh.af.Tnil/Papers/pdfyae9l.pdf 


5a=n,  +  n2H - +«„_].  (3.14) 

The  bitwise  and  operation  is  denoted  here  by  the  symbol  A. 
The  symbol  =*  denotes  a  bitwise  barrel  roll  to  the  right.  That 
is,  “s=>j”  means  shift  the  integer  s  to  the  right  by  j  digits. 
Hence,  the  result  of  the  operation  “(s/\2a)=>a”  is  either  1 
or  0,  depending  on  whether  or  not  a  particle  occupies  the  ath 
local  state.  Notice  that  the  symbols  a *  and  a  are  overloaded, 
so  that  when  they  are  used  with  a  single  argument,  that  ar¬ 
gument  is  interpreted  as  a  spin  value.  If  a*  and  a  are  used 
with  two  arguments,  the  first  argument  is  interpreted  as  a 
spin-index  and  the  second  argument  is  interpreted  as  a  ket. 

Notice  that  these  symbolic  definitions  of  the  multiple-spin 
creation  and  annihilation  operators  use  the  basis-state  symbol 
T'  on  the  LHS  of  the  rules,  but  "9  is  not  used  at  all  on  the 
RHS  in  the  definition  of  the  rules.  Hence,  it  may  seem  that 
the  use  of  the  symbol  is  superfluous  here.  However,  this  is 
not  the  case,  because  its  use  allows  me  to  define  the  action  of 
the  creation  and  annihilation  operators  on  a  superposition  of 
basis  states  in  a  recursive  fashion: 

a't[a,A"i#[s'\  +  B~\= A  <A[a',1P[,s,]]+at[a:,.6],  (3.15) 

a[a>A'¥[s'\  +  3]=A  a[o',*F[.s,]]+ a[a,B],  (3.16) 

Using  this  convention,  it  is  possible,  for  example,  to  destroy 
a  spin  in  local  state  a  of  a  superposed  state,  say  T'f.Sj] 
by  directly  supplying  this  state  as  the  second  argu¬ 
ment.  Then,  Eq.  (3.13)  correctly  expands  out  to 
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a[a,'^[5j]  +  ^[j2]]=a[a,'?'[5'i]]+a[a,^[52]]. 

(3.17) 

If  the  special  symbol  1Jf  were  not  used,  then  one  would  get 
the  wrong  answer, 

a[a,sl+s1\=a[a,${\,  (3.18) 

where  s3  is  the  numerical  sum  of  and  s2.  Of  course,  it  is 
possible  to  use  a  special  symbol  in  place  of  the  plus  sign  to 
represent  superimposed  states.  I  have  chosen  not  to  do  this. 
With  the  symbol  convention,  mathematica  can  by  de¬ 
fault  manipulate  expressions  involving  the  superposition  of 
an  arbitrary  number  of  states  and  represent  them  in  memory 
in  a  compact  fashion.  After  the  action  of  the  collision  opera¬ 
tor  (which  is  mathematically  defined  earlier  in  this  paper  and 
symbolically  defined  immediately  below)  on  to  a  superposed 
state,  the  resulting  new  state  in  general  has  identical  basis 
states  that  are  repeated  in  the  superposition,  where  each  oc¬ 
currence  has  a  different  amplitude.  Using  the  ^  symbol  con¬ 
vention,  all  these  types  of  replications  are  automatically  re¬ 
duced  down  to  the  one  term,  since  mathematica 
automatically  adds  coefficients  of  common  terms. 

Next,  the  multiple-spin  number  operator  is  defined  as  a 
composition  of  the  multiple-spin  creation  and  annihilation 
rules 

nfa,’Ir]=a^t«,a[a,''If]].  (3.19) 

With  rules  (3.12),  (3.13),  and  (3.19),  for  the  creation,  anni¬ 
hilation,  and  number  operators,  it  is  then  straightforward  to 
implement  the  universal  gate,  Eq.  (3.3),  by  composition: 

U[si  ,s2,^r]  =  '^-Cnt[s2  ,a[j!  ,,P']]-itat[si  ,n[^2^]] 

+  (A  - 1  )«|>ia  ,¥]  +  (D-  1  )«[s2  ,¥] 
-(AJrD)n[sun[s2,V\\  (3.20) 

where  the  c  numbers  A,  B,  C,  and  D  are  components  of  an 
SU(2)  matrix  (£  2). 

In  the  case  of  the  quantum  lD3Px  model,  the  collision 
operator  mixes  the  on-site  kets,  [Oil)  and  [  1 00) .  Three 
qubits  are  affected.  I  use  a  modified  rule  to  directly  handle 
this  situation.  The  on-site  collision  operator  for  the  lD3Px 
quantum  lattice  gas  is  implemented  by  the  following  compo¬ 
sition  of  universal  gates: 

>*ib  ,S2>^]=¥-  Cnf[s2  ,a[sla  ,a[slft , •'?]]] 

-Ba\su,a\slb,a[s2^Vi\ 

-(l-/f)n[jla,n[s16,-^]] 

-(i-£>)»[x2,n 

+  (1-D)«[ji0,«[j2,^]] 

+  (1-£>)«[52,«[5i*,^]] 

-(A -D)n[sla  ,n[s16  ,n[s2  ,¥]]]. 

(3.21) 
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The  lattice-gas  collision  operator  according  to  Eq.  (2.7) 
for  the  V—lf  system  is  thus  defined  as  a  sevenfold  compo¬ 
sition 

C[^]=I/t20,21,19,U[17,18,16,I/[14,15,13,U[ll,12,10, 

X  U[  8,9,7,C/[5,6,4,  U[  2,3,1  ]]]]]]] .  (3.22) 

This  is  actually  handled  recursively  in  the  symbolic  imple¬ 
mentation,  so  C  works  regardless  of  the  size  of  the  system. 

The  streaming  operator  for  the  quantum  lattice  gas  is 
implemented  using  two  rules,  one  to  stream  the  right  moving 
particles,  denoted  S+  ,  and  the  other  to  stream  the  left  mov¬ 
ing  particles,  denoted  .  Note  that  the  right  moving  par¬ 
ticles  occupy  local  states  2, 5, 8,1 1,14,17,20  and  the  left  mov¬ 
ing  particles  occupy  local  states  3,6,9,12,15,18,21.  5+  and 
S_  are  defined  in  terms  of  a  sevenfold  composition  of  inter¬ 
change  operators  (3.5): 

5+m=Y[2,5,Y[5,8,x[8,ll,Y[ll,14,Y[14,17, 

X*[  17,20,*]]]]]],  (3.23) 

5-m=Y[21,18,Y[18,15,Y[15,12,^[12,9,Y[9,6, 

x*[6, 3, *]]]]]].  (3.24) 

Again,  these  are  handled  recursively  in  the  symbolic  imple¬ 
mentation,  so  the  streaming  operators  work  regardless  of  the 
size  of  the  system.  A  global  shift  of  particles  is  done  by 
successive  local  interchanges  of  particles  occupancies  [18], 
Finally,  the  evolution  rule,  denoted  E,  for  the  entire  quan¬ 
tum  system  is  the  composition  of  the  last  three  rules 

£[*]=C[S+[S..[*]]].  (3.25) 

Any  other  compound  rules  that  may  be  needed  in  a  simula¬ 
tion  can  be  defined  in  a  similar  fashion  by  composing  pre¬ 
defined  simpler  rules.  Therefore,  beginning  with  a  superposi- 
ton  of  basis  states  'b(/)  =  2t<?!>?*[.?]  the  dynamical  evolution 
equation  corresponding  to  Eq.  (2.3)  is 

<h(r+r)  =  E[‘l>(t)],  (3.26) 

where  the  result  is  a  new  superposition  over  a  different  set  of 
basis  states  <J>(f +t)=2s/(£s#*[.s']. 

B.  The  lD3Px  model 
1.  Classical  version 

Let  us  consider  a  simple  lattice-gas  model  as  a  concrete 
example,  called  the  lD3Px  lattice-gas  model,  in  this  paper. 
This  model  was  first  studied  by  Qian  in  1990  [19]  and  is 
referred  to  as  Model  I  in  his  thesis.  The  lattice  gas  is  one 
dimensional  and  has  three  bits  per  site,  a  rest  particle  with 
mass  two  and  speed  ±  1  particles  with  mass  one.  The  mass 
and  momentum  at  a  lattice  site  is 

m  =  2nQ+nx+n2  and  px=nx-n2.  (3.27) 

There  are  two  local  configurations  both  with  m= 2  and  px 
=  0:  (1)  {n0,«i,«2}={lAO}  and  (2)  {n0,nx  ,«2}={0,1,1}. 
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m  =2 


head-on 

— 


rest 


FIG.  1.  Head-on  collision  in  the  lD3Px  lattice-gas  model.  The 
single  equivalence  class  has  m  =  2  and  px=  0. 


Kfl(x  +  /ea><+r)  =  na(x,t)H-fia(«>l<).  (3.31) 

For  the  lD3Px  model,  the  lattice  vectors  are  e0 =0,  e,  =x, 

and  e2 — ~x  and  the  collision  term  is  specified  by  the  single 
function 


These  two  configurations  are  members  of  the  only  collision 
set  (which  is  called  an  equivalence  class).  An  equivalence 
class  has  two  or  more  members.  Figure  1  illustrates  the 
equivalence  class  of  the  lD3Px  model.  Its  two  elements  are 
the  configuration  of  two  head-on  particles  {01 1}  and  the  con¬ 
figuration  with  a  single  rest  particle  {100}. 

Because  the  total  number  of  particles  and  the  total  mo¬ 
mentum  must  be  conserved,  the  collision  part  of  the  dynam¬ 
ics  can  only  permute  the  local  configurations.  The  collision 
equation,  which  is  applied  homogeneously  across  the  lattice, 
can  be  expressed  as  in  terms  of  a  mapping  function  U  as 
follows: 


s'  =  U(s),  (3.28) 

where  U  maps  2s  configurations  to  2s  new  configurations. 
For  the  simple  lD3Px  lattice,  U  is 

U{{  011})={100}, 

U({100})  =  {011}. 

If  a  configuration  s  is  not  a  member  of  an  equivalence  class, 
then  U(s)  =s.  In  other  words,  if  the  incoming  state  is  not  a 
member  of  an  equivalence  class,  then  the  outgoing  state  is 
set  equal  to  the  incoming  state.  To  speed  up  a  lattice-gas 
simulation,  the  mapping  function  U  may  be  precomputed 
before  the  simulation  and  accessed  in  lookup  table  fashion 
during  the  simulation. 

In  a  computer  implementation,  it  is  convenient  to  use  two 
arrays  to  simultaneously  store  the  states  s  and  s'.  Therefore, 
in  Eq.  (3.28),  data  in  file  array  that  stores  the  “incoming” 
state  s  is  transformed  by  the  action  of  the  lookup  table  U 
(which  is  applied  homogeneously  over  the  entire  array)  and 
the  output  is  written  into  the  next  array  to  store  the  new 
“outgoing”  state  s' . 

It  is  conventional  to  write  the  collision  rule  in  terms  of  the 
occupation  variables  na  —  1  or  0,  which  are  Boolean  values. 
The  collision  rule,  expressed  for  an  individual  local  state,  is 
written 


n'a(x,t)=na(xyt)+n,a{n *),  (3.29) 

where  the  collision  term  fM«*)  =  ±l  or  0.  Writing 
ila(n*)  with  an  asterisk  subscript  on  a*  denotes  that  the 
collision  term  for  the  uth  local  state  depends  on  all  the  on¬ 
site  local  states.  It  is  conventional  to  write  the  streaming  rule 
in  terms  of  na  also, 

na(x+/ea,t+r)=n'a(x,t).  (3.30) 

Combining  Eqs.  (3.29)  and  (3.30),  the  microscopic  transport 
equation  is  therefore 


n=/ii«2(l  -/20)-n0(l  -»,)(1  — n2),  (3.32) 

where  O0=fi  andfilj2=  -fi.  Then  explicitly  for  the  lD3Px 
model,  the  microscopic  transport  equation  (3.31)  is 

n0(x,t+r)  =  n0(x,t)  +  fl(x,t),  (3.33) 


nia(x±/>t+T)  =  nh2(x,t)~-Cl(x,t). 

A  lattice-Boltzmann  equation  describes  the  dynamics  of  the 
lD3Px  lattice-gas  system  at  the  mesoscopic  scale.  The  me¬ 
soscopic  average  of  the  occupation  variable  na(x,t)  is  the 
probability  of  occupancy 

fa(x,t)^(na(x,t)).  (3.34) 

Here,  the  angle  brackets  around  a  microscopic  quantity  de¬ 
note  its  mesoscopic  expectation  value  obtained  by  ensemble 
averaging.  The  kinetic  transport  equations  are 

fo(x,t+T)=f0(x,t)  +  {n(x,t)),  (3.35) 

/u(x±/,t+T)=/u(x,0-(f2(x,0). 

To  carry  out  a  classical  lattice-gas  simulation  at  the  mesos¬ 
copic  scale,  we  can  approximate  fimESO(x,/)s(fl(x,0)  by  a 
mean-field  collision  term,  denoted  f2mf(x,t),  that  neglect 
particle-particle  correlations: 

<ft(x,/)>~n^x,0«/i/2(l  -/0)-/o(l  ~/i)(l  -/2). 

(3.36) 

A  statement  of  detailed  balance  can  be  written  by  setting  the 
mean-field  value  of  the  collision  term  (3.36)  to  zero  at  equi¬ 
librium 


(H)-nmf(O=0.  (3.37) 

Therefore,  the  probability  of  occupancies  satisfies  the  equa¬ 
tion 


foq= 


/r/?+d  -/TOO  -JT)‘ 


(3.38) 


This  equation,  along  with  equations  for  the  mass  and  mo¬ 
mentum  densities 


Po=2fi0q+fiq+fV  and  ux0=f?-fq,  (3.39) 

gives  us  a  nonlinear  system  of  three  equations  in  five  un¬ 
knowns/^,  po,  and  ux0  ■  Hence,  it  is  possible  to 

analytically  solve  for  the  occupation  probabilities  flq,  f)q, 
and /I*1  in  terms  of  p0  and  ux0 .  When  the  system  is  at  rest 
at  equilibrium,  px=0,  then  fiq=f2q=d  and  the  probability 
of  occupancy  for  the  rest  particle  state  is 


046702-9 


JEFFREY  YEPEZ 


PHYSICAL  REVIEW  E  63  046702 


f?= 


l-2d+2d 2 


(3.40) 


Using  Eq.  (3.36),  the  Jacobian  of  the  collision  Jab 
is 


/ 


J= 


- 1  +  2d-2d2 


\-2d+2d2 


(1  ~d)d 

(1  ~d)d 

1  -2d+2d2 

1  —2d+2d2 

(d-  l)d 

(d-l)d 

1  —2d+2d2 

1  -2d+2d2 

(d-l)d 

(d-l)d 

1  —2d+2d2 

1  -2d+2d2 

The  eigenvectors  of  J  are 

1 1 ) — (2,1,1 ), 
|2)  =  (0,1,—  1), 
<(l-2d+2d2)2 


|3>  = 


did- 1) 


,1,1  • 


(3.41) 

(3.42) 

(3.43) 

(3.44) 


The  eigenvectors  |l)  and  |2),  corresponding  to  mass  and 
momentum,  span  a  two-dimensional  hydrodynamic  sub¬ 
space.  The  remaining  eigenvector  |3)  is  a  kinetic  eigenvec¬ 
tor,  which  in  this  case  is  density  dependent.  The  eigenvalues 
of  J  are 


o 

II 

(3.45) 

>" 

NJ 

l! 

o 

(3.46) 

1  —  2<7+6<72  — 8<73  +  4<74 

(3.47) 

- 1  +2d-2d2 

Now  using  the  lattice  vectors  e0=0,  ex-l,  and  e2 
=  - 1,  and  the  expression  for  J  given  in  Eq.  (3.41),  we  set 
the  secular  determinant  of  the  linearized  Boltzmann  equation 
equal  to  zero 

[(et(/,'*°'k+aT)-l)8ab-Joh]=0.  (3.48) 

This  allows  us  to  solve  for  the  dispersion  relations  for  the 
lattice-gas  system  obeying  what  is  called  generalized  hydro¬ 
dynamics.  Equation  (3.48)  is  a  result  from  the  generalized 
hydrodynamics  of  classical  lattice-gas  systems  previously 
worked  out  by  Das,  Bussemaker,  and  Ernst  [20]  and  Grosfils, 
Boon,  Brito,  and  Ernst  [21].  Taking  /=  r=  1 ,  we  get  the 
following  dispersion  relation: 

(l-2d+2d2)e3a-2[d-3d2+4d3-2d4 

+  (1  —3d+3d2)cosk]e2w  (3.49) 

+  (l-2d)2ll+2d(d-l)cosk]ea+4d2(d-\)2=0. 

(3.50) 


0  1  2  3  4  5  6 

Wave  Number  k  ( 2n/IS) 


Wave  Number  k  (2 K/ty 

FIG.  2.  The  real  part  of  the  dispersion  relation  for  the  mesos¬ 
copic  lD3Px  lattice  gas  in  the  long  wavelength  limit  and  mean-field 
limit  at  a  reduced  background  density  of  <7=0.214  286. 

This  is  a  cubic  equation  in  e“,  and  it  is  analytically  solvable. 
The  only  hydrodynamic  mode  is  a  damped  sound  wave 
a>(k)=±csk+iT(p)k2.  Real  and  imaginary  parts  of  the  dis¬ 
persion  relations  for  the  lD3Px  lattice-gas  model  are  shown, 
respectively  in  Fig.  2  and  Fig.  3.  The  real  part  of  the  disper¬ 
sion  relations  indicates  a  sound  mode  [Re(rn)— >±csk  as  k 
— >0].  The  imaginary  part  of  the  dispersion  relation  for  the 
hydrodynamic  mode  is  parabolic  for  small  wave  numbers, 
indicating  viscous  damping  of  the  sound  mode  [Im(<w) 
-*Tk2  as  k—>0].  The  sound  damping  constant  T  approaches 
zero  as  the  background  mass  density  approaches  zero  [19]. 
That  is,  low-mass  density  waves  can  oscillate  without  vis¬ 
cous  damping. 

The  real  part  of  the  dispersion  relation  for  the  sound  mode 
for  the  lD3Px  lattice-gas  model  set  with  a  background  den¬ 
sity  of  </=6/4K,  with  F=7/,  is  shown  in  Fig.  2.  The  real 
part  of  the  dispersion  relation  indicates  a  sound  mode 
[Re(a>)— >±csk  as  £— > 0  where  c5=0.74//r].  The  data 
points,  plotted  as  black  circles,  are  solutions  to  the  linearized 
Boltzmann  equation  in  the  mean-field  limit.  The  curves  with 
slope  of  ±  cs  are  numerical  linear  fits  to  the  data.  The  imagi¬ 
nary  part  of  the  dispersion  relation  for  the  sound  mode  for 


Wave  Number  k 

FIG.  3.  The  imaginaiy  part  of  the  dispersion  relation  for  the 
mesoscopic  lD3Px  lattice  gas  in  the  long  wavelength  limit  and 
mean-field  limit  at  a  reduced  background  density  of  <7=0.214286. 
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the  lD3Px  lattice-gas  model  is  shown  in  Fig.  3.  The  imagi¬ 
nary  part  of  the  dispersion  relation  indicates  sound  damping 
[Im(o))— >fT&2  as  £-+0  where  T  =  0.08 The  parabola 
is  a  numerical  fit  to  the  data  in  the  region  of  small  k<  1 .  The 
calculations  shown  in  Figs.  2  and  3  were  done  with  a  mass 
density  filling  fraction  of  do=6/4F=0.214,  where  a  small 
system  size  of  K=7/  is  used.  In  this  case,  k=2'ir/V 
=  0.898. 

2.  Quantum  version 

A  hypothetical  lattice-based  quantum  computer  (with 
computational  sites  depicted  as  circles)  arranged  as  a  one¬ 
dimensional  lattice  is  shown  in  Fig.  4.  At  each  lattice  site 
resides  B= 3  qubits  in  ID  in  this  example  with  V—lf  sites. 
The  on-site  ket  \  if/)  resides  in  a  25-dimensional  submanifold. 
The  large  circle  on  the  right  represents  an  expanded  view  of 
this  on-site  submanifold,  which  is  denoted  by  B.  The  basis 
states  of  B  are  shown  in  the  number  representation.  Each  site 
is  coupled  to  its  nearest  neighboring  sites  by  a  mechanism 
allowing  for  the  exchange  of  qubits.  If  the  exchange  mecha¬ 
nism  retains  all  quantum  entanglement  (and  thereby  spread¬ 
ing  it  through  the  quantum  computer),  then  the  quantum 
computer  is  considered  fully  coherent.  If  the  exchange 
mechanism  is  classical  (destroying  quantum  entanglement  by 
collapsing  the  wave  function),  then  it  is  called  a  type  II  quan¬ 
tum  computer  (which  is  simply  a  large  array  of  small  quan¬ 
tum  computers  interconnected  by  a  classical  communication 
network). 

The  associated  lD3Px  quantum  lattice-gas  model  has 
three  qubits  per  site,  |<7a)  =  «fl|0)  +  /?Jl)  for  a  =  0,l,2.  The 
zeroth  qubit  represents  a  rest  particle  of  mass  two  and  the 
first  and  second  qubits  represent  moving  particles  of  speeds 
±  1 ,  translating  in  the  right  and  left  going  directions,  respec¬ 
tively. 

The  m-  2,  px= 0  equivalence  class  is  spanned  by  the 
states  1 100)  and  |01 1).  Collisional  entanglement  occurs  only 
between  these  two  states,  £|100)  +  *|011),  where  £  and  x 
are  c  numbers.  The  on-site  ket,  \iff)-\qv)®\qi)®\qi),  is 

l^)=AAA|Hl>+)8ol9ia2|110>  +  ^o«iA|10I> 

+/8o«.«2|100)  +  «0AA|011>  +  aoA«2|010) 

+  a0a1y32|001)  +  a0a1«2|000).  (3.51) 

The  outgoing  on-site  ket  \if/')  =  U\tf/)  is 


fioPiPi 
PoPl<*2 

aPoaia2  +  ba0p1p2 
ca0BiB2  +  da0l3i/32 
a0A«2 

«O<*102 
«0“l«2 

1  0  0  0  0  0  0  0\  IPoPiPi 

0  1  0  0  0  0  0  01  P<sPx(*i 

0  0  1  0  0  0  0  0  PbOLyp-2 

OOOabOOO  /?o“i«2 
OOOcdOOO  ctQpip2  ’ 

0  0  0  0  0  1  0  0  aQpla2 

0  0  0  0  0  0  1  0  a0«iA> 

0000000  1/  a0atia2l 

(3.52) 

where  the  local  collision  operator  is  the  8X8  matrix  with 
one  2X2  block,  which  is  a  member  of  the  U(2)  unitary 
group  satisfying 

|a|2+|h|2=|c|2+|d|2=l,  (3.53) 

ac*  +  bd*=a*c+b*d= 0,  (3.54) 

|a|2  +  |c|2  =  |b|2  +  |</|2=l,  (3.55) 

ab*+cd*=a*b  +  c*d= 0.  (3.56) 

The  quantum  lD3Px  lattice  gas  obeys  detailed  balance 
because  the  collision  operator  0  is  a  unitary  matrix  [8]. 

The  mass  and  momentum  densities  for  the  quantum 
lattice-gas  system  are 

P=2{qQ\n\q0)+{qx\n\ql)-k{q2\n\q2),  (3.57) 

Ux={d\\n\qi)-(q2\n\q2).  (3.58) 

Viscous  dissipation  does  not  necessarily  occur  in  quantum 
lattice-gas  systems.  Global  entanglement  of  the  wave  func- 
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FIG.  5.  Damping  of  a  mass  density  wave  for  a  system  with  V—l  sites  in  the  classical  lD3Px  model  simulated  using  a  mesoscopic 
Boltzmann  equation  with  the  collision  term  expressed  in  the  mean-field  approximation.  The  background  density  is  rfo  =  6/4E=0.214.  The 
ordinate  is  the  absolute  value  of  the  amplitude  of  the  mass-density  wave  divided  by  the  peak  amplitude  of  the  initial  perturbation. 


tion  significantly  complicates  the  dispersion  relations,  which 
are  determined  by  the  following  equation: 


Det 


x)Sah 


d{ye%U'naU-na)®\®  ■  ■  -011^) 
— 


(3.59) 


where  1Peq  is  the  steady-state  equilibrium  wave  funtion, 
which  is  the  ground  state  of  the  system.  I  have  explicitly 
written  the  collision  operator,  as  in  Eq.  (2.23),  in  spatially 
separated  form.  In  general,  as  described  in  Sec.  II  Bl,  a 
—  aQ  +  a,  where  aQ  is  an  index  that  refers  to  the  first  local 
state  at  some  particular  site  in  the  system.  According  to  the 
ordered  numbering  scheme  used,  «o=0  at  the  first  site  of 
the  system,  ao~B  at  the  second  site,  a0=2B  at  the  next 
site,  and  so  on.  Without  loss  of  generality,  in  Eq.  (3.59)  we 
can  assume  we  are  working  at  the  first  site  of  the  system 
where  na=na®  1®  •  •  •  ®  1.  In  the  classical  case,  C  is  a  per¬ 
mutation  matrix  and  the  steady-state  equilibrium  wave  func¬ 
tion  is  a  tensor  product  over  the  on-site  kets 

v 

|^eq)=  ®  |^).  (3.60) 

X=1 


In  turn,  the  on-site  kets  are  formed  by  a  tensor  product  over 
the  individual  qubits 


l<T>=  ®  (^lO+VWTIO)).  (3.61) 

a~\ 


Finally, and/^q=ai2/[d2  +  (l  ~d)2]  according  to  Eq. 
(3.40).  The  Jacobian  of  the  collision  matrix  element  appear¬ 
ing  in  Eq.  (3.59)  is  computable  using  Eqs.  (3.60)  and  (3.61) 
[see  Eq.  (3.41)  in  Sec.  Ill  B  1],  In  the  quantum  mechanical 
case,  |1Teq)  is  not  expressible  as  a  tensor  product  state,  and 
hence  the  Jacobian  of  the  collision  matrix  element  appearing 
in  Eq.  (3.59)  becomes  complicated. 


C.  Simulations 
1.  Classical  simulation 

A  time  history  of  the  mass  density  wave  for  a  small  sys¬ 
tem  with  F=7/  sites  is  shown  in  Fig.  5.  The  exponential 
envelope  is  analytically  determined  by  an  analysis  of  the 
linearized  lattice-Boltzmann  equation  in  the  mean-field  limit 
(see  Fig.  3).  The  predicted  sound  damping  constant  T 
=  0.08/2/r  is  in  excellent  agreement  with  the  simulation 
data. 

Plotted  in  Fig.  6  are  damping  time  constants  of  mass  den¬ 
sity  waves  in  the  classical  lD3Px  lattice  gas  for  different 
system  sizes  from  V—2 /  up  to  V-256/.  The  log-log  plot 
shows  the  power-law  behavior,  known  as  diffusive  ordering , 
typical  of  a  lattice-gas  system  in  the  viscous  regime.  The 
power  law  in  this  case  is  7’=0.44F2,  which  is  parabolic. 
Each  circle  is  determined  from  a  mesoscopic  scale  simula¬ 
tion  that  was  initialized  with  a  sinusoidal  perturbation  of 
8p=0.04m/f  from  a  uniform  background  mass  density  at 
half-filling,  p=2/w//.  The  damping  constant  T =/2/r  is  de¬ 
termined  from  the  envelope  of  the  resulting  standing  wave 
e~llTcos  cot  (see  Fig.  5).  The  mean-field  estimates  of  the 
damping  time  constant  are  the  circles.  The  line  is  a  linear 
best  fit  to  these  estimates.  The  estimated  damping  constant 
deviates  only  slightly  from  power-law  behavior  at  the  small¬ 
est  system  sizes.  This  is  an  example  of  “fluidlike”  behavior 
occurring  in  systems  far  below  the  continuum  limit.  The  in¬ 
set  plot  is  a  linear  plot  of  the  data  for  F=£  16  and  the  parabola 
is  the  same  diffusive-ordering  power-law  in  the  larger  log- 
log  plot. 


2.  Quantum  simulation 

I  have  tested  the  quantum  lattice-gas  formalism  described 
in  this  paper  by  carrying  out  exact  numerical  simulations  of 


FIG.  6.  Diffusive  ordering  in  the  classical  lD3Px  model  com¬ 
puting  at  the  mesoscopic  scale  using  the  mean-field  approximation. 
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FIG.  7.  Mass  and  momentum  sectors  of  the  lD3Px  lattice-gas 
model  with  V=l/  sites  plotted  versus  the  number  of  states  per 
sector. 

a  lD3Px  model,  which  is  described  in  detail  in  Sec.  Ill  B  1. 
In  this  section,  I  present  results  obtained  from  the  numerical 
simulation  of  a  small  system  with  V=lf  sites.  I  have  used 
the  symbolic  numerical  technique  described  above  in  Sec. 
Ill  A  2.  The  principal  finding  is  that  the  quantum  lattice  gas 
does  not  display  viscous  damping. 

Since  the  evolution  operator  conserves  mass  and  momen¬ 
tum,  we  can  divide  the  Hilbert  space  into  disjoint  mass- 
momentum  sectors.  When  the  lattice-gas  evolution  operator 
maps  a  particular  state  residing  in  a  mass-momentum  sector 
to  a  new  state,  the  new  state  must  also  reside  in  that  same 
mass-momentum  sector.  The  Hilbert  space  for  the  F=7/ 
system  has  over  two  million  dimensions.  The  number  of 
states  within  each  mass-momentum  sector  of  the  V=lf  sys¬ 
tem  are  graphically  illustrated  in  Fig.  7.  The  density  plot  on 
the  left  side  of  Fig.  7  clearly  shows  that  the  allowable  mass- 
momentum  sectors  are  all  contained  within  a  hexagonal 
boundary.  The  distribution  for  the  number  of  available  states 
within  a  mass-momentum  sector  is  reflection  symmetric 
about  half-filling  (m  — 14)  and  about  zero  momentum  (pr 
=  0). 

I  have  simulated  the  F=7/  system  (with  BV=  21  glo¬ 
bally  entangled  qubits)  in  the  mass  m  =  6  and  momentum 
Px~ 0  sector.  In  this  mass-momentum  sector,  there  are  5376 
basis  states.  The  goal  of  the  numerical  test  was  to  measure 
the  sound  damping  constant  in  the  quantum  lD3Px  model 
and  compare  the  result  to  the  mean-field  estimate.  The  sys¬ 
tem  was  initialized  with  a  sinusoidal  perturbation  of  the  mass 
density  field,  with  a  wavelength  equaling  the  grid  size  of  the 
periodic  system  (A  =  F).  All  the  states  in  the  m= 6,  px= 0 
sector  were  superposed  by  choosing  amplitudes  in  such  a 
fashion  that  the  entropy  of  the  initial  state  is  maximized, 
subject  to  the  independent  constraints  of  conservation  of 
probability,  mass,  and  momentum.  The  entropy  function  was 
taken  to  be 

S=~J,  [|ctt|2ln|cj2+(l-|cj2)ln(l-|cj2)], 

(3.62) 

where  ca  is  the  amplitude  of  the  ket  |a)  in  the  m  =  6,  px 
=  0  mass-momentum  sector.  Given  a  particular  desired  pro¬ 
file  of  the  mass  density  field,  it  is  more  difficult  to  construct 
an  initial  state  that  completely  resides  in  only  one  sector  than 
tb  use  an  initial  state  that  spans  the  entire  Hilbert  space. 


Site  ® 


FIG.  8.  Initial  mass  density  sinusoidal  perturbation  in  the  quan¬ 
tum  lD3Px  lattice  gas  for  a  small  system  size  of  V-  7/  with  peri¬ 
odic  boundary  conditions.  The  total  number  of  qubits  in  the  simu¬ 
lation  is  BV=  21.  The  simulation  is  initialized  with  a  sinusoidal 
perturbation  in  the  m  —  6,  px=0  mass-momentum  sector  with  a 
peak  amplitude  of  Sp— 0.4  from  a  uniform  background  mass  den¬ 
sity  at  po=  2g  =  0.857.  So  the  fractional  mass  density  variation  is 
initially  one  part  in  two,  which  is  an  extremely  large-scale  fluctua¬ 
tion.  The  wavelength  equals  the  system  size.  The  initial  mass  den¬ 
sity  field  is  not  exactly  sinusoidal,  because  aside  from  the  limitation 
of  only  V=7/  sites,  it  is  produced  by  the  interference  of  all  5376 
in  the  m  =  6  and  px=  0  sector.  An  algorithm  using  Lagrangian  mul¬ 
tipliers  maximizes  the  entropy  of  the  resulting  wave  function  and 
chooses  all  the  amplitudes  of  the  initial  state. 


However,  it  is  computationally  advantageous  to  limit  the 
simulation  to  a  single  sector  of  the  Hilbert  space,  so  that 
memory  allocation  in  the  computer  is  kept  at  a  manageable 
level.  Figure  8  shows  a  maximized  entropy  state  used  in  the 
test  simulation  presented  in  this  section. 

The  data  from  the  simulation  run  is  presented  in  several 
ways.  First,  the  peak  amplitude  of  the  mass  density  wave  is 
recorded  after  every  time  step.  The  amplitude  is  normalized 
in  such  a  fashion  that  at  time  f=0  it  has  unity  value.  In  a 
viscous  fluid  with  sound  damping,  the  peak  amplitude  would 
oscillate  and  decay  exponentially  in  time  by  the  factor, 

e  „  cos(2tc,j///),  where  cs  is  the  sound  speed  and  T  is  a 
positive  definite  damping  constant  as  is  shown  in  Fig.  5. 
However,  for  the  quantum  lD3Px  model,  the  numerical  re¬ 
sult  indicates  F  may  be  zero  for  certain  collision  operators. 

A  time  series  history  of  the  square  of  the  peak  amplitude 
is  plotted  in  Fig.  9,  using  the  same  format  as  Fig.  5  for  the 
classical  lD3Px  model  with  the  same  grid  size  and  initial 
condition.  In  the  quantum  simulation,  the  peak  amplitude 
does  not  decay  in  time,  unlike  the  results  obtained  in  the 
classical  lattice-Boltzmann  simulations  shown  in  Fig.  5.  Ini¬ 
tially,  within  the  first  couple  of  dozen  time  steps,  the  peak 
amplitude  appears  to  decay,  very  much  like  it  does  in  a  clas¬ 
sical  microscopic  simulation  or  lattice-Boltzmann  simulation 
of  the  lD3Px  model.  However,  the  amplitude  does  not  con¬ 
tinue  to  damp  in  subsequent  time  steps.  The  peak  amplitude 
rises  and  falls.  No  damping  is  observed  even  after  a  thousand 
time  steps.  An  expanded  view  of  the  first  250  time  steps  is 
shown  underneath.  Since  the  algorithm  is  unitaiy  (and  hence 
the  collisions  obey  the  principle  of  detailed  balance)  the  dy¬ 
namics  is  reversible. 

In  Fig.  10,  these  data  are  presented  in  scatter  plot  fashion, 
where  the  square  of  the  normalized  peak  amplitude  is  plotted 
versus  its  first  order  time  derivative.  I  used  the  following 
difference  formula  to  approximate  the  time  derivative: 
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FIG.  9.  Oscillations  of  a  mass  density  wave  in  the  quantum 
lD3Px  lattice  gas  for  a  system  size  of  K=7/  in  the  m—6  and  px 
—  0  sector.  The  ordinate  is  the  absolute  value  of  the  amplitude  of 
the  mass-density  wave  divided  by  the  peak  amplitude  of  the  initial 
perturbation. 

Sp\x,t)  p2(x,t—  t)  —  p2(x,t+  t) 

The  data  appear  randomly  scattered,  but  is  clustered  along  a 
“cone”  corresponding  to  the  speed  of  sound  in  the  lD3Px 
model,  which  the  Boltzmann  analysis  of  Sec.  Ill  B  1  predicts 
to  be  c,=0.74//r. 

To  obtain  a  more  accurate  estimate  of  the  sound  speed  in 
the  quantum  lD3Px  simulation,  a  Fourier  transform  of  the 
time  series  history  of  the  mass  density  at  a  single  site  of  the 
system  was  computed  and  the  power  spectrum  p*(x)plo(x) 
plotted  (see  the  bottom  plot  of  Fig.  1 1).  The  top  plot  shows 
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FIG.  10.  Normalized  peak  (absolute  value  of  the  amplitude  of 
the  mass-density  wave  divided  by  the  peak  amplitude  of  the  initial 
perturbation)  versus  the  first  derivative  of  the  normalized  peak  of 
oscillations  of  a  mass  density  wave  in  the  quantum  lD3Px  lattice 
gas  for  a  system  size  of  V—l/  in  the  m  =  6  and  px=  0  sector.  We 
have  plotted  maximum  speed  curves  corresponding  the  individual 
particle  velocity,  c=±/ It.  As  expected,  all  the  data  are  contained 
within  this  “cone.”  In  addition,  we  have  plotted  sound-speed 
curves  corresponding  toct=±0.74//r,  which  is  analytically  deter¬ 
mined  from  a  mean-field  approximation  of  the  system  using  the 
linearized  lattice-Boltzmann  equation.  Most  of  the  data  is  clustered 
around  the  sound-speed  curves,  and  additional  data  points  scattered 
within  the  “sound-speed  cone”  indicates  randomness  in  the  oscil¬ 
lation  of  the  mass  density  wave. 
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FIG.  11.  Time  history  of  the  mass  density  at  site  x=6f  for  a 
system  with  V=1  /  sites  plotted  versus  time.  A  discrete  Fourier 
transform  of  this  time  series  data  is  taken  to  give  p*(x)pwfr).  A 
peak  in  the  power  spectrum  |pj2  occurs  at  about  0.72/7 r,  which  is 
close  to  the  expected  sound  speed.  The  abscissa  is  converted  into 
unit  of  velocity,  c=/7r,  to  show  that  there  is  a  unique  sound  speed. 
The  ordinate  has  units  of  (m/fx  t)2. 

the  time  series  collected  by  measuring  the  fluctuation  of  the 
mass  density  field  of  the  F=7/  quantum  lD3Px  lattice-gas 
system.  The  signal,  which  is  p(6/,t),  is  measured  at  site  x 
=6/.  Plotted  below  is  the  power  spectrum  of  the  Fourier 
transform  of  the  signal,  which  is  |pj2,  versus  sound  speed 
(this  is  proportional  to  the  oscillation  ffequeny,  cs-/f).  A 
peak  in  the  power  spectrum  occurs  just  below  the  mean-field 
approximation  of  sound  speed,  cs  -  0.1  A/ I  r,  which  is  plotted 
as  the  vertical  bar.  (See  Fig.  2  for  the  mean-field  value  esti¬ 
mate  of  cs.) 

IV.  CONCLUSION 

The  main  results  of  this  paper  are  as  follows. 

The  quantum  mechanical  wave  equation  is  recast  as  a 
lattice-Boltzmann  equation  describing  a  quantum  lattice-gas 
system. 

The  continuity  and  Navier-Stokes  equations  constitute  a 
macroscopic  effective  field  theory  for  the  quantum  lattice- 
gas  system  and  quantum  entanglement  changes  the  value  of 
the  transport  coefficients. 

A  symbolic  math  method  was  presented  for  simulating 
dynamical  quantum  systems. 

With  reversible  microscopic-scale  dynamics,  a  feature  of 
classical  lattices  is  that  dissipation  occurs  at  the  macroscopic 
scale.  However,  viscous  damping  is  not  observed  in  simula¬ 
tions  of  the  quantum  lD3Px  lattice-gas  model,  which  is  also 
microscopically  reversible. 

The  sound  speed  of  mass  density  waves  is  the  same  as  the 
classical  value. 

Given  the  memory  and  speed  constraints  of  classical  com¬ 
puters,  today  only  small  quantum  lattice  gas  can  be  exactly 
simulated.  I  have  performed  many  test  simulations  of  the 
quantum  lD3Px  model  for  system  sizes  ranging  from  V 
=3/  up  to  V—l/  and  have  included  results  from  the  V 
—  1/  quantum  simulation  in  the  paper,  since  this  was  the 
largest  computer  run. 

I  do  not  wish  to  argue  that  results  obtained  for  such  a 
small  system,  with  V=l/  sites,  can  give  us  too  much  in¬ 
sight  about  the  true  macroscopic  behavior  of  the  quantum 
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lattice  gas,  which  is  only  well  defined  in  the  continuum  limit. 
Further  testing  is  required  on  larger  systems  and  in  two  and 
three  dimensions  and  will  be  presented  in  a  subsequent  pa¬ 
per.  Yet,  in  the  classical  version  of  the  model,  hydrodynami¬ 
clike  behavior  is  observed  in  very  small  systems  (see  Figs.  5 
and  6).  The  type  of  behavior  found  in  the  small  V= 7/  quan¬ 
tum  lattice-gas  system  may  also  occur  for  larger  systems.  So, 
quantum  lattice  gases  of  multiple  grid  sizes  should  be  simu¬ 
lated.  To  this  end,  a  compiled  version  of  quantum  lattice-gas 
code  is  being  developed  in  fortran  90  and  will  be  tun  on 
available  supercomputers. 

The  issue  of  the  similarity  or  distinction  between  particle- 
particle  correlations  (as  occurs  in  classical  lattice  gases)  and 
quantum  entanglement  (as  occurs  in  quantum  lattice  gases) 
has  not  been  addressed  in  this  paper.  Yet,  this  is  an  issue  that 
can  be  studied  using  quantum  lattice-gas  simulations. 
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APPENDIX  A:  DERIVATION  OF  THE  QUANTUM 
LATTICE-GAS  TRANSPORT  EQUATION 
IN  THE  CONTINUUM  LIMIT 

In  this  appendix,  I  would  like  to  rederive  the  transport 
equation  (2.19)  for  the  quantum  lattice-gas  system.  The  deri¬ 
vation  given  here  is  carried  out  in  the  continuum  limit  (imag¬ 
ine  a  space-time  lattice  with  infinite  resolution  as  the  cell 
size  vanishes).  All  the  usual  restrictions  arising  from  the  dis¬ 
cretization  of  the  microscopic  quantities  are  temporarily  re¬ 
moved.  A  particle  can  exist  at  any  point  in  space  and  time, 
and  it  can  also  have  any  momentum  p  =  mv.  The  only  as¬ 
sumption  I  make  here  is  that  I  can  still  decompose  the  space 
time  into  an  ordered  set  of  local  states,  which  in  this  case  is 
infinite  but  denumerable.  That  is,  I  imagine  there  are  an  in¬ 
finite  number  of  local  states  at  each  point  in  space  (2?=°°), 
one  corresponding  to  every  possible  particle  momentum. 
Since  the  number  of  points  in  the  space  is  also  infinite  ( V 
=  <»),  the  total  number  of  local  states  are  doubly  infinite 
(N=BV=«>2).  Nevertheless,  I  assume  the  local  states  are 
well  ordered  and  denumerable. 

The  probability  of  finding  a  particle  with  momentum  p  in 


the  ath  local  state  located  at  position  x  given  by  Eq.  (2.1 1)  is 
the  following  matrix  element: 

/(xV,0^(^(t)l«J^(t)).  (Al) 

I  assume  f(x,p,t)  is  a  continuous  and  differentiable  mesos¬ 
copic  field  quantity.  For  the  moment,  suppose  the  a  is  the 
local  state  of  an  “incoming”  particle,  preceding  a  possible 
collision  event.  I  still  want  to  imagine  the  particle  dynamics 
divided  into  mutually  exclusive  events  (collision  followed  by 
streaming)  repeated  in  stepwise  fashion  ad  infinitum.  Next, 
the  probability  of  finding  a  particle  in  the  local  state  a', 
corresponding  to  momentum  p'  at  position  x'=x 
+(r/m)p',  is  expressed  by  the  matrix  element 

f(x  +  (r/m)p'  ,p'  ,t+T)^(^(t+T)\na,\'i'(t+T)). 

(A2) 

Suppose  a'  is  the  local  state  of  the  “outgoing”  particle. 
Then,  a  basic  definition  of  the  total  time  derivative  of 
f(x,p,t)  is  the  following  ratio: 

df(x,p,t)  f(x  +  (r/m)p'  ,p'  ,t)-f(x,p,t) 


or,  in  terms  of  the  matrix  elements,  it  is 

df(x,p,t ) 
dt 

W* +  T)|n„,|*('+  r)) -<¥(0|« J¥(*)> 

=  lim - . 

T — *0  7 

(A4) 

This  is  the  seed  of  a  Boltzmann  equation  for  particle  trans¬ 
port  and  the  RHS  of  this  equation  constitutes  the  collision 
term,  although  this  may  not  appear  quite  obvious  at  this  point 
in  the  development.  In  the  following  development,  I  shall 
interpret  the  collision  term  and  rewrite  it  so  that  it  explicitly 
depends  only  on  na  at  positions  and  I'FQ)).  In  so  doing,  we 
shall  see  how  the  collision  dynamics  is  inherently  encoded  in 
this  expression. 

First,  we  add  zero  to  the  RHS  of  the  above  equation  to 
write  the  collision  term  in  two  parts,  explicitly  separating  the 
total  change  into  “temporal-change”  and  “spatial-change” 
parts,  as  follows: 

df{x,p,t) 

dt 

,.  Wf  +  r)|n„,|^(f+r))-(^(f)|«a,|^)) 

=  lim - 

T — >0  7 

+  lim - . 

T— >0  T 

(A5) 
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From  the  time-displacement  operation,  eT3l3‘f(x,p,t) 
= f(x,p,t+  r),  we  see  that  the  first  term  on  the  RHS  of  the 
above  equation  is  a  partial  derivative  with  respect  to  time 

df(x+{T/m)p,p,t )  , 

— — +  0(Sh2) 

at 

{^{t+r)\na,\^{t+T))-{nt)\na'\nt)) 

-  hm - . 

T— +0  T 

(A6) 

The  Stouhal  number,  Sh,  is  defined  as  the  ratio  of  the  mean- 
free  time  to  the  characteristic  length  scale  (Sh=  r/t).  Simi¬ 
larly,  from  the  space-displacement  operation,  en,  vf(x,p,t) 
= f(x+ w,p,t),  we  see  that  the  second  term  is  a  partial  de¬ 
rivative  with  respect  to  position 

v-  Vf(x,p,t)  +-(v- V)2f(x,p,t)  +  0(Kn3) 

_  ]Im - 

T— >0  ^ 

(A7) 

The  Knudsen  number,  Kn,  is  defined  as  the  ratio  of  the 
mean-free  path  to  the  characteristic  length  scale  (Kn 
~/IL).  Therefore,  we  have  the  convective  derivative 

df(x,p,t )  df(x+(rm/p),p,t )  ,  ,  . 

— Tt - - T, - 

+  j  ($‘  v)V(jc»P»0  +  0(Sh\Kn3), 

(A8) 


This  result  is  expected,  since  in  quantum  mechanics,  the  par¬ 
tial  time  derivative  of  an  operator  is  found  by  calculating  the 
commutator  of  that  operator  with  the  Hamiltonian.  Using  this 
result,  the  Boltzmann  equation  (A5)  becomes 


df(x,p,t)  _  i 
dt  h 


I  nat  —  n„  \ 

+  lim  ¥(*)  — - ^  ty(t)  .  (A1 1) 

T— >0  \  T  / 


Now  the  RHS  no  longer  depends  on  [^(f  +  t))  (so  it  is  local 
in  time),  but  it  is  still  nonlocal  in  space  because  it  depends  on 
na<  as  well.  That  is,  if  the  RHS  of  the  above  equation  were 
to  depend  only  on  na,  then  it  would  have  “strictly  local” 
form. 

Third,  using  the  fact  that  e‘llrlfl=SC,  we  can  rewrite  the 
commutator  as 


Now,  na  and  na<  are  related  by  the  similarity  transformation 
(2.14),  so  the  commutator  reduces  to 


-\_nai  ,//]=  lim 


&naC-na, 


Inserting  this  into  Eq.  (A1 1)  gives  the  final  local  form  of  the 
quantum  Boltzmann  equation  for  f(x,p,t),  which  is 


composed  of  a  local  term  and  a  nonlocal  advection  term.  In 
the  local  term,  it  is  technically  correct  (albeit  unconven¬ 
tional)  to  explicitly  write  the  partial  time  derivative’s  depen¬ 
dence  on  r,  even  though  r— >0.  This  is  done  to  stress  an 
equivalence  with  the  matrix  element  formulation  given  by 
Eq.  (A5). 

Second,  we  rewrite  the  “local  change”  term.  Since 
|^(r+T))=e‘"T/ft|^(0>  and  eikrlfl=  1  +itir/h  +  0(i2), 
we  have 

('^(t  +  T)|na.|1?(t+T)) 

«<¥(/)|jv|¥(/)>+ jmt)\[na<  Mnt)) 

+  0{t1).  (A9) 

Using  this  equation  in  conjuction  with  Eq.  (A6),  we  have 

dflx+^p,p,t) 

h~A~Jt - (A10) 


df(x,p,t)  1  .  „ 

--  /  -  «Hm-(TO|C?«flC-i.Jg(f)).  (A14) 

T— *0  ' 

Notice  that  the  collision  term  depends  only  on  the  wave 
function  evaluated  at  time  f  and  the  occupancy  of  the  ath 
local  state  located  at  position  x.  However,  if  there  exists 
quantum  superposition  between  particles  at  different  points 
in  space,  then  |^(0)  cannot  be  written  in  separable  tensor 
product  form  over  the  spatial  points.  So  in  this  case,  the 
collision  term  is  “nonlocal.”  Hence,  when  I  say  the  lattice- 
Boltzmann  equation  is  local  in  form,  I  mean  this  in  a  pseudo- 
classical  sense,  barring  nonlocal  quantum  entanglements: 
And  this  is  why  I  said  in  the  introduction  of  this  paper  that 
the  lattice-Boltzmann  equation,  which  accounts  for  global 
entanglement  through  the  collision  process,  is  an  exact  refor¬ 
mulation  of  the  many-body  Schrodinger  equation. 

There  is  one  more  point  to  make  in  this  appendix.  From 
the  basic  definition  (A3)  for  the  total  time  rate  of  change  of 
f(x,p,t),  we  see  that  Eq.  (A14)  can  be  written  as  the  follow¬ 
ing  “finite-difference”  equation 
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f(x+(T/m)p',p',t)=f(x,p,t)  +  ('¥(t)\ClnaC-na\V(t)). 

(A15) 

This  is  the  lattice-Boltzmann  equation  [see  Eq.  (2.19)  in  Sec. 
II B  2].  It  is  important  to  note  that  the  Boltzmann  equation  is 
still  an  exact  representation  of  the  particle  dynamics,  even 
when  expressed  in  finite-difference  form.  This  is  immedi¬ 
ately  obvious  when  the  identity  na=S^na,S  is  inserted  into 
the  collision  term,  ^(t)^  S^nat§d-na\ty(t)).  Then,  the 
lattice-Boltzmann  equation  becomes  a  simple  identity 


f(x+{Tlm)p'  ,p'  ,t)=f{x,p,t)  +  (^{t+  T)\na,\^{t+  t)) 
-(¥(0l«J*«)>.  (A16) 

In  the  case  of  a  finite  resolution  lattice  (used  in  a  computa¬ 
tional  simulation  of  the  quantum  lattice-gas  system),  the 
lattice-Boltzmann  equation  is  the  appropriate  formulation  of 
the  particle  dynamics.  However,  the  quantum  Boltzmann 
equation  (A14),  in  differentiable  point  form,  becomes  the 
appropriate  formulation  of  the  particle  dynamics  when  talk¬ 
ing  about  the  system  in  the  continuum  limit. 


APPENDIX  B:  REPRESENTATION  OF  A  TWO-QUBIT 
GATE  FOR  A  2-SPIN  SYSTEM 

In  this  appendix,  we  show  that  Eq.  (3.3)  is  a  manifestly 
unitary  operator  that  entangles  two  qubits  according  to  the 
SU(2)  special  unitary  group.  Let  us  consider  a  quantum  spin 
system  with  only  two  spins.  Then  the  Hilbert  space  is  four 
dimensional,  and  we  choose  the  following  basis  kets  in  the 
number  representation: 


ll\ 

l°\ 

l°\ 

l°\ 

0 

1 

0 

0 

0 

,  |10>  = 

0 

,  |oi) = 

1 

.  lll>= 

0 

\o] 

\oj 

\0/ 

\i/ 

(Bl) 


In  this  basis,  the  creation  operators  are 


0 

0 

0 

1° 

0 

0 

°\ 

-1 

0 

0 

0 

>  a2  = 

'  0 

0 

0 

0 

0 

0 

0 

°l 

i 1 

0 

0 

0 

0 

0 

-1 

oj 

\0 

1 

0 

0/ 

(B2) 

Since  a\  and  a\  have  real  components,  the  annihilation  op¬ 
erators  are  the  transpose  of  the  matrices  given  in  Eq.  (B2): 
a i ~(a\)T  and  al={a\)T.  The  universal  gate  operator  is  ex¬ 
pressed  in  terms  of  the  following  five  operators: 


a\a2  = 


and 


1° 

0 

0 

°\ 

1° 

0 

0 

°\ 

0 

0 

-1 

0 

f  0 

0 

0 

0 

0 

0 

0 

0 

,  a\ai  = 

-1 

0 

0 

u 

0 

0 

oj 

u 

0 

0 

0 

«l(l~"2)  = 


(l-«i)«2  = 


/  0  0  0  0\ 
0  10  0 
0  0  0  0 
(0  0  0  0/ 

/  0  0  0  0\ 
0  0  0  0 
0  0  10 
\o  0  0  0 / 


1  — «1  —  «2  = 


/  10  0 
0  0  0 
0  0  0 
0  0  0 


0  \ 
0 
0 


1  / 


(B3) 


(B4) 


We  can  represent  a  block  diagonal  4X4  unitary  matrix  in 
terms  of  these  five  operators  as  follows: 


I  1  0  0 

0  A  B 
0  CD 
\0  0  0 


0 

0 

0 

-1/ 


—Anl(l-n2)-Ba\a2-Calal 


+  D(l-n1)n2+l-nl-n2. 

(B5) 


When  the  2  X  2  block  is  a  member  of  SU(2)  as  given  by  Eq. 
(3.1),  this  expression  for  a  unitary  matrix  becomes  a  repre¬ 
sentation  of  a  universal  gate  given  by  Eq.  (3.3). 

In  this  appendix,  we  used  a  two-spin  quantum  system  as 
an  example  system  for  illustrating  how  a  universal  gate  can 
be  expressed  in  terms  of  the  multispin  creation  and  annihila¬ 
tion  operators.  Although  we  used  a  two-spin  system  in  this 
example,  the  procedure  outlined  here  also  works  for  a  spin 
system  with  an  arbitrary  number  of  spins. 

All  permutations  of  single  fermion  states  may  be  imple¬ 
mented  by  successive  application  of  a  “interchange  opera¬ 
tor”  [22],  here  denoted  by  xap'  >  where  the  permutations 
occur  between  state  a  at  site  x  and  states  /?'  at  site  x' 


Xap’ => a\a p, +(1^010+1- ala e-a^, a p, .  (B6) 


n  1 


*+  '"'t 

—  n  1  n  —  n  1 


This  is  a  special  case  of  the  universal  quantum  gate,  Y 
where  0=  77-/2,  £=0  and  £=0.  The  interchange  operator 
correctly  handles  any  necessary  phase  change  due  to  the  anti¬ 
commutation  relations  (3.2). 
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